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ABSTRACT 

A  numerical  model  which  employes  observed  stratospheric  winds  to 
advect  simulated  tracers  was  developed.  This  model  successfully  repro- 
duces many  qualitative  features  of  the  observed  fields  of  both  ozone 
and  radioactive  debris.  As  the  tracers  evolve,  the  horizontal  eddies 
constitute  the  principal  process  modifying  the  tracer  zonal  mean. 
Although  the  vertical  eddies  and  the  zonal  mean  cell  are  both  an  order 
of  magnitude  weaker  than  this  process,  the  latter  of  these  tends  to 
always  act  in  the  same  sense  so  that  its  effect  becomes  more  important 
over  long  periods  of  time. 
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I.  BACKGROUND 

A.  OBSERVATIONS 

The  study  of  trace  substance  transport  in  the  stratosphere  can 
provide  many  useful  insights  into  the  dynamics  of  that  region.  The 
history  of  an  evolving  inert  tracer  depends  upon  its  initial  state  and 
upon  several  statistical  properties  of  the  wind  field.  In  the  strato- 
sphere, detailed  observations  of  the  wind  field  are  not  always  plentiful 
enough  to  reveal  these  rather  subtle  characteristics  of  the  general 
circulation.  In  addition  to  observations  of  the  wind,  there  exist  some 
observations  of  the  distributions  of  such  natural  tracers  as  ozone  and 
radioactive  debris.  These  observations  may  be  used  to  extend  the  avail- 
able data  base,  and  theories  of  the  general  circulation  may  be  tested 
by  requiring  them  to  simulate  the  advection  of  natural  tracers. 

The  influence  of  stratospheric  transport  processes  was  first  noticed 
when  the  early  spectrophotometer  observations  of  Dobson  and  his  co-workers 
(Dobson  et  al .  1928)  revealed  that  the  distribution  of  total  atmospheric 
ozone  has  the  following  characteristics: 

1.  The  total  amount  of  ozone  per  unit  area  increases  from  the  equator 
toward  the  poles. 

2.  In  high  latitudes  this  quantity  varies  with  season  exhibiting  a 
pronounced  maximum  in  late  winter  or  early  spring. 

3.  No  apparent  correlation  exists  between  total  ozone  and  solar 
activity. 

4.  Total  ozone  is  correlated  with  the  passage  of  upper-level  synoptic 
systems,  attaining  positive  departures  from  its  long-term  mean  in  the  lows 
and  negative  departures  in  the  highs. 
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Dobson's  observations  are  consistent  with  those  obtained  by  the  81 
station  network  in  the  Soviet  Union  and  reported  by  Bozkov  (1968).  These 
observations  indicate  that: 

1.  Minimum  total  ozone  occurs  at  10  N. 

2.  Maximum  total  ozone  occurs  at  60  to  70  N,  with  values  decreasing 
toward  both  the  equator  and  the  pole. 

3.  The  maximum  zonal  mean  gradient  of  total  ozone  occurs  between  30 
and  60  N. 

4.  The  total  ozone  tends  to  have  a  wavelike  distribution  in  the 
horizontal  with  low  values  occurring  over  the  western  portions  of  conti- 
nents and  high  values  over  the  eastern  portions. 

The  photochemical  theory  of  ozone  formation  (See  Craig,  1965)  pre- 
dicts that,  at  photochemical  equilibrium,  the  distribution  of  total 
ozone  should  attain  a  maximum  in  low  latitudes  at  the  subsolar  point 
and  increase  toward  the  poles.  Brewer  (1949)  and  Dobson  (1956)  attempt 
to  reconcile  this  theory  with  observation  by  postulating  a  mean  cell  with 
ascending  motion  in  low  latitudes,  northward  advection  in  the  high  mid- 
latitude  stratosphere,  descent  near  the  pole,  and  a  return  circulation 
at  low  levels  in  mid-latitudes.  This  circularion  supposedly  carries  the 
ozone,  which  is  generated  by  solar  ultraviolet  radiation  in  low  latitudes, 
northward  and  downward  into  the  lower  polar  stratosphere  where  it  ac- 
cumulates since  there  it  is  protected  from  photodissociation. 

The  Brewer-Dobson  theory  is  seriously  challenged  by  airborne  filter 
observation  of  the  distribution  of  Tungsten-185  from  the  Hardtack  nuclear 
tests  (Feely  and  Spar,  1960).  The  center  of  concentration  of  this  debris 
remained  at  50  mb  and  nearly  at  the  latitude  of  injection.  The  cloud  did 
not  migrate  northward  and  downward,  but  rather  spread  in  the  meriodinal  plane 
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in  such  a  manner  that  the  isopleths  of  concentration  sloped  down- 
ward from  the  equator  toward  the  pole. 

B.  PARAMETERIZED  MODELS 

Reed  (1950)  and  Normand  (1953)  note  that  stratospheric  troughs  are 
regions  of  both  high  temperature  and  high  total  ozone.  This  leads  them 
to  postulate  that,  since  both  potential  temperature  and  ozone  mixing 
ratio  increase  upward,  the  lows  must  be  regions  of  descending  motion 
and  the  highs  regions  of  ascent.  Furthermore,  the  troughs  and  ridges 
move  more  slowly  than  the  wind  so  that  horizontal  advection  shifts  the 
centers  of  greatest  mixing-ratio  departure  downwind  from  the  associated 
pressure  system.  Because  the  wind  almost  always  has  a  westerly  component 
in  the  winter  season,  the  ozone  maxima  lie  in  the  southerly  winds  east 
of  the  troughs  and  the  minima  are  east  of  the  ridges  in  northerly  winds. 
Thus,  the  departures  of  both  temperature  and  mixing-ratio  are  explained 
in  terms  of  the  systematic  interaction  between  horizontal  and  vertical 
advection.  High  mixing-ratio  and  temperature  is  accompanied  by  descent 
and  poleward  motion  while  low  values  of  these  two  quantities  occur  with 
rising  and  equatorward  motion. 

Citing  earlier  work  by  Dodson  (1960),  Newell  (1961)  carries  this  line 
of  reasoning  one  step  further.  He  explains  the  poleward  slope  of  the 
isopleths  by  invoking  eddy  mixing  in  which  northward  motion  is  correlated 
with  descent  and  southward  motion  with  ascent.  In  the  same  paper  he 
computes  the  flux  of  ozone  by  correlating  total  ozone  with  the  winds  at 
12-18  km  in  the  maximum  ozone  layer.  He  finds  that  90  percent  of  his 
computed  flux  is  due  to  large-scale  eddies  and  that  the  total  flux  is 
more  than  adequate  to  explain  observed  changes  in  the  distribution. 
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Prabhakara  (1961)  generates  zonal  mean  ozone  distributions  which 
depend  upon  the  interaction  of  photochemistry,  a  mean  cell  that  descends 
near  the  pole,  and  large-scale  eddies  parameterized  as  anisotropic  dif- 
fusion. His  results,  which  are  in  good  qualitative  agreement  with  obser- 
vation,show  that  such  transport  processes  can  explain  many  departures 
from  photochemical  equilibrium. 

In  contrast  to  Prabhakara' s  work,  in  which  the  principal  axis  of 
eddy  diffusion  is  horizontal,  there  have  been  a  number  of  models  simu- 
lating the  transport  of  radioactive  debris  with  anisotropic  diffusion 
schemes  in  which  the  principal  axis  is  inclined  to  the  horizontal  (Reed 
and  German,  1965;  Davidson,  et  al . ,  1965;  Seitz,  et  al . ,  1968;  and  Fair- 
hall  and  Reed,  1968).  In  such  models  it  is  assumed  that  the  northward- 
downward  velocity  correlation  in  the  large-scale  eddies  can  be  para- 
meterized in  terms  of  a  diffusion  tensor  whose  principal  axis  slopes 
downward  toward  the  pole.  This  means  that  the  tracer  is  presented  with 
a  path  of  least  resistance  parallel  to  the  observed  slope  of  the  iso- 
pleths  of  mixing  ratio  in  the  real  atmosphere.  In  order  to  keep  the 
tracer  on  these  sloping  paths  and  to  prevent  excessive  vertical  spreading, 
it  is  necessary  to  make  the  vertical  component  of  the  diffusion  tensor 

much  smaller  than  that  along  the  principal  axis.  Davidson  et  al .  (1966) 

9   2    -1 
use  a  principal  diffusivity  of  4  x  10  cm  sec   and  a  vertical  dif- 

4  2    -1 
fusivity  of  about  10  cm  sec  .  However,  since  characteristic  horizontal 

3 
distances  are  on  the  order  of  10  km  while  vertical  distances  are  about 

1  km  and,  since  the  speed  of  diffusion  depends  upon  the  square  of  these 

distances,  a  ratio  of  10  between  the  two  diffusivities  is  not  physically 

unreasonable. 
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C.  GENERAL  CIRCULATION  MODELS 

Hunt  and  Manabe  (1968)  and  Hunt  (1969)  report  the  most  extensive  and 
realistic  simulation  of  stratospheric  tracers  to  date.  In  their  experi- 
ments various  simulated  substances  are  introduced  into  the  wind  fields 
generated  by  a  general  circulation  model  and  the  history  of  the  numeri- 
cally evolved  mixing  ratio  fields  are  examined. 

In  the  1968  experiments  inert  tracers  are  permitted  to  evolve  from 
zonal ly  symmetric  initial  states.  The  first  such  tracer  roughly  corre- 
sponds to  radioactive  debris  from  an  equatorial  injection  and  is  initially 
distributed  in  a  band  extending  from  the  equator  to  ten  degrees  north  with 
the  maximum  concentration  at  50  mb  pressure  height. 

The  second  experiment's  initial  state  represents  the  distribution  of 
ozone  at  photochemical  equilibrium  throughout  the  region  of  integration. 
Both  tracers  are  treated  as  inert  substances  inasmuch  as,  except  for 
removal  of  any  tracer  that  finds  its  way  into  the  troposphere,  no  sources 
or  sinks  are  included  in  these  experiments. 

The  simulated  histories  of  both  tracers  are  reassuringly  similar  to 
observations  of  real  tracers  in  the  atmosphere.  Both  the  radioactive 
debris  and  the  ozone  begin  at  once  to  migrate  northward  and,  to  some 
extent,  downward.  Thus  the  concentration  at  the  equator  decreases  and 
that  at  the  poles  increases,  so  that  the  initial  gradient  becomes  smaller 
and  eventually  reverses. 

The  actual  transfer  by  the  quasi-horizontal  eddies  is  well  illustrated 
by  the  evolution  of  the  mixing  ratio  field  at  various  horizontal  levels  in 
the  model.  The  zonally  symmetric  character  of  the  initial  field  quickly 
disappears  becoming  sinusoidal  with  wave  number  four  predominating  as  it 
does  in  the  circulation  generated  by  the  model.  With  the  passage  of  time 
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the  waves  becomes  increasingly  exaggerated,  forming  the  mixing-ratio 
contours  into  long  northward  protrusions.  Eventually  the  protrusions 
separate,  leaving  isolated  islands  of  high  concentration  in  northern 
latitudes. 

Hunt  (1969)  reports  two  additional  experiments  with  the  same  model 
and  involving  ozone  transport  with  photochemical  sources  and  sinks  added 
The  mechanism  of  transfer  in  these  studies  and  in  those  previously  dis- 
cussed are  very   similar.  In  both,  the  mean  cell,  which  produces  a 
convergence  of  tracer  into  the  subtropics,  is  opposed  by  an  eddy  diver- 
gence out  of  that  region,  with  the  eddies  being  particularly  effective 
in  transferring  tracer  northward  in  the  region  poleward  of  30N.  The 
principal  difference  introduced  by  photochemistry  is  the  strong  source 
in  low  latitudes  which  prevents  reduction  or  reversal  of  the  initial 
poleward  gradient  by  replacing  the  ozone  as  fast  as  transport  processes 
can  remove  it. 

The  weak  link  in  any  experiment  involving  either  parameterization 
or  a  general  circulation  model  lies  in  the  arbitrary  specification  of 
the  transport  process  and  its,  at  best,  indirect  relation  to  the  actual 
motions  in  the  stratosphere.  In  view  of  these  difficulties,  it  seems 
instructive  to  compare  advection  of  simulated  tracers  by  actual  observed 
winds  with  both  observation  and  previous  attempts  to  model  both  ozone 
and  radioactive  debris.  The  present  study,  therefore,  employs  actual 
observed  stratospheric  winds  to  advect  simulated  tracers  in  the  presence 
of  sub-grid-scale  diffusion  and  numerically  integrates  the  equation  of 
continuity  for  the  tracer  in  time  to  obtain  the  evolution  of  an  inert 
tracer  from  an  arbitrary  initial  state. 
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This  approach  is  largely  independent  of  any  preconceived  ideas  about 
the  stratospheric  general  circulation  and  so  provides  a  third  source  of 
information  intermediate  between  the  parameterized  models  or  the  general 
circulation  models  and  observations  of  real  tracers. 

II.  MODEL 

A.  BASIC  EQUATION 

The  present  study  follows  the  lead  of  Hunt  and  Manabe  (1968)  by 
numerically  exploring  the  detailed  advection  of  trace  substances  with- 
out resort  to  arbitrarily  defined  diffusion  tensors.  Furthermore,  the 
winds  used  to  accomplish  this  advection  are  based  on  actual  observations 
of  the  real  atmosphere  and  are  not  the  result  of  a  mathematical  model 
that  may  or  may  not  reproduce  the  properties  of  the  actual  general 
circulation. 

The  basic  equation  on  which  this  study  depends  is  the  equation  of 
continuity  for  an  inert  tracer  without  sources  or  sinks,  expressed  in 
flux  form  and  spherical  pressure  coordinates. 

where  m  is  the  mixing  ratio  of  the  tracer,  R  the  earth's  radius,  0  the 
geographic  latitude,  7\   the  longitude,  p  the  pressure,  QO  the  substantial 
pressure  derivative  following  an  individual  air  parcel,  and  u  and  v  the 
wind's  eastward  and  northward  components  respectively. 

The  local  mean  of  any  dependent  variable  in  (1)  is  defined  as  its 
average  computed  over  a  volume  characteristic  of  the  grid  distance  to  be 
used  in  the  integration.  The  value  at  any  point  within  the  region  can 
now  be  represented  by  the  sum  of  this  mean  and  a  departure  from  the  mean 
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(2) 

V  -  V  +  V  *  co-  CO  +  cO   ^ 

where  the  starred  variables  are  the  departures  and  the  barred  ones  are 
the  means.  These  relations  are  then  substituted  into  (1)  and  the  entire 
relation  is  averaged  recognizing  that,  while  the  average  of  a  fluctuation 
alone  or  of  the  product  of  a  fluctuation  with  a  mean  vanish,  the  averaged 
product  of  two  fluctuations  does  not.  (See  Haltiner  &  Martin,  1957). 

The  terms  containing  products  of  fluctuations  are  then  approximated 
by  using  the  Prandtl  mixing-length  assumption  that  the  correlation  between 
velocity  and  mixing-ratio  is  proportional  to  the  mean  gradient  of  mixing 
ratio  as  given  by  __ 


-i^^ .  (4) 
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The  vertical  and  horizontal  eddy  diffusivities ,  K..  and  K?  are  assumed 
to  be  constant  in  both  space  and  time.  This  particular  choice  of  con- 
stants, with  the  same  eddy  diffusivity  in  both  horizontal  directions  and 
without  any  off-diagonal  elements  in  the  diffusivity  tensor,  implies  that 
diffusion  is  isotropic  in  the  horizontal  plane  with  the  axis  of  least 
diffusion  oriented  vertically  since  K,«K?. 

By  substituting  (4)  for  the  eddy  correlation  terms  in  (3)  and  writing 
the  flux  terms  back  in  advective  form  through  use  of  the  mass  continuity 
equation  for  the  local  mean  flow,  one  obtains 


'&CD  -  co  '$5£>    Ar    —   QCQ    +     ca         3  no 

dp*-       f£-cas>fid<pK       T  &PJ      /Lco$*<ft     S>         (5) 


At  this  point  (5)   is  placed  in  non-dimensional ized  form  by  the  fol- 
lowing linear  transformation  of  the  velocity  fields: 

&=Ctot{'j     \7^V0V/>        g5  -   CJ>0co'  .  (6) 

Here  the  sub-zero  values  are  constants  with  magnitudes   (u     =  v     =20  kt, 

oo 
_i 
Co)    =  5  mb  day     )  so  chosen  that  the  primed  variables,  whose  variations 

contain  the  spatial   and  temporal   changes  of  the  velocities,  are  of  order 

one.     Since  the  equation  is  linear  in  m  no  such  scaling  needs  to  be  done 

for  thatvariable  and  (5)  becomes 

-*  ^m      ^f--1 -*£&   ,     '      ^rr}?       (7) 

For  solution  on  a  digital  computer  the  infinitesimal  derivativies  in 
(7)  have  to  be  replaced  by  finite  differences,  thus  converting  the  partial 
differential  equation  into  a  system  of  algebraic  equations  in  the  mixing 
ratio  and  wind  components  at  discrete  points  in  the  region  of  integration. 
This  procedure  requires  that  the  independent  variables  of  the  problem  be 
expressed  as  multiples  of  small  but  finite  increments  so  that  the  co- 
ordinates of  a  point  (p,$,?l,  t)  become  i,j,k  and  n?  where  i  =  p/£p, 
j  =  {i>  -  &)/&<#  ,  k  =^/S/\  and  n  =  t/  £>t  and  the  values  assumed  by 
(p>0»^»  t)  are  required  to  make  i,  J,  k  and  n  integers.  In  this  for- 
malism m(p,0 ,9^ ,  t)  ism"  ,  .  and  v(p,^,^,  t)  is  v"  .  .. 


Once  such  a  coordinate  system  is  set  up  it  becomes  possible  to 
approximate  the  time  derivative  by  a  forward  difference 

3t  St        J  (8) 

the  first  order  space  derivatives  by  centered  differences,  for  example, 

sp  2?^  Jip  >         (9) 

and  to  use  the  three  point  approximation  in  place  of  the  second  order 
space  derivatives 

The  A  operator  always  applies  at  the  point  i,  j,  k;  although  the  space 
subscripts  are  usually  omitted  to  simplify  notation  motion.  The  sub- 
script on  the  operator  denotes  the  coordinate  direction  along  which  the 
difference  is  to  be  taken  and  the  superscript  the  order  of  that  dif- 
ference. 

With  all  the  derivatives  replaced  by  finite  differences  and  after 
multiplication  through  by&t,  (7)  becomes 


where  the  n  and  n  +  1  superscripts  denote  the  base  time  level  and  the 
newly  generated  time  level  respectively,  and  the  $   superscript,  which 
also  denotes  time,  can  take  on  only  one  or  the  other  of  these  two  values 
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Each  of  the  factors  in  brackets  in  (14)  is  the  product  of  St 
with  a  collection  of  constants  which,  when  combined,  has  the  units  of 
reciprocal  time.  Thus  each  bracket  forms  the  dimensionless  ratio  of 
the  timestep  to  a  time  characteristic  of  the  particular  process  chang- 
ing the  mixing  ratio.  For  example,  in  the  first  term  2£p/ca;  represents 

2 
the  vertical  advection  time  or  in  the  fourth  term  (^p)  /K-.  is  the 

meridional  diffusion  time. 

Since  each  of  these  factors  is  a  constant,  (11)  is  simplified  by 

replacing  them  with  single  quantities  to  reproduce  the  simpler  form 

(12) 


The  pressure  increment  is  12.5  mb,  that  of  latitude  5°  and  that  of 
longitude  10°.  With  this  information  the  first  three  advection  times 
are  evaluated  to  give: 

vs=  2izs^/u0  =  2.5"  days.  (13) 

The  value  of  VL  does  not  represent  the  actual  advection  time  enter- 
ing into  the  calculation  since  a  factor  of  cos/  is  not  included  in  its 
evaluation.  This  factor  represents  the  decrease  in  advection  time  as 
the  distance  corresponding  to  6?v  shrinks  with  approach  to  the  pole. 
Since  the  cosine  of  the  latitude  varies  from  .819  at  35  N  to  .087  at 
85  N  the  actual  zonal  advection  time  changes  by  nearly  an  order  of 
magnitude  over  the  region  of  integration. 
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Since  a  state  of  nearly  total  ignorance  exists  with  regard  to  the 
actual  magnitudes  of  the  eddy  diffusivities  in  the  stratosphere, the 
diffusion  time  along  a  given  axis  is  set  to  a  simple  multiple  of  the 
advection  time  along  the  same  axis.  In  this  study  the  multiple  (called 
l/b7)  is  chosen  to  be  eight.  This  means  that  the  parameterized  dif- 
fusion can  remove  departures  from  equilibrium  1/8  as  fast  as  advection 
operating  by  itself.  This  assumption  leads  to  the  following  values  for 
the  diffusion  times: 

%  =  s»o  days  .  04) 

The  difference  of  a  factor  of  four  between  TV  and  "Q  arises  because 

'kTX   =  2S^and,  since  the  horizontal  diffusion  is  assumed  to  be  isotropic, 

the  transport  must  occur  four  times  as  rapidly  over  half  the  distance. 

The  values  of  K,  and  K?  corresponding  to  T.   and  T).  are  1  x  10  cm 
sec"  and  4.5  x  10  cm  sec"  compared  to  2  x  10  and  1  x  10   used  in 
Prabhakara's  parameterized  study.  That  author's  definition  of  diffusion 
and  the  one  used  here  are  essentially  different.  In  the  former  case  the 
diffusion  was  used  to  simulate  transfer  by  the  large-scale  eddies,  while 
in  this  study  it  has  only  to  account  for  transfer  by  fluctuations  too 
small  to  be  represented  in  synoptic-scale  data.  The  similarity  between 
the  two  vertical  diffusivities  is  unfortunate;  one  would  have  preferred 
to  use  a  vertical  diffusivity  that  differed  from  that  of  Prabhakara  by 
a  factor  of  at  least  ten  because,  in  the  parameterized  studies,  diffusion 
is  called  upon  to  provide  much  of  the  transport  that  advection  accomplishes 
in  this  model . 

The  entire  treatment  of  diffusion  is  rather  artificial  but  since  it 
serves  as  a  computational  smoother,  its  effect  can  not  be  minimized  with- 
out sacrificing  numerical  stability. 
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B.  TIME  INTEGRATION 

Integration  in  time  is  accomplished  using  a  backward  corrected 

Euler  scheme  (see  Kurlhara,  1965),  If  the  advective  and  diffusive 

terms  in  (12)  may  be  replaced  by  f(m  ,  ,  v11,  un)  this  scheme  may  be 

represented  thus: 

rr/-  nV  +  f  OT)°.  co"  V'7  <x°  )  4 

J    >      ;  (15) 

;     y       *  (16) 

Given  the  field  of  m  and  the  wind  at  the  n  time  level  m*-  is  computed 
from  (15).  Then  m^  and  the  winds  at  the  n  +  1  time  level  are  used  in 
(16)  to  obtain  m   .  Thus  m    is  first  predicted  using  a  forward 
difference  and  then  corrected  with  a  backward  difference,  the  entire 
process  being  repeated  at  each  successive  time  level  until  m  is  known 
for  the  entire  period  of  interest. 

According  to  Kurlhara  (1965),  if  this  method  is  applied  to  the  one 
dimensional  advection  equation, 

(18) 

stability  will  result  for  b-6'1.  By  experiment,  it  was  discovered  that 

optimum  results  are  obtained  by  using  a  time  step  £>t  ■  1/48  day  producing 

the  following  values  for  the  b's  (b7  =  1/8) 

b,  =  4.i7*io3,  b^  =  i.  67*10*   b3-  £ 33*lo  3 

b««  0.55*15^  65-4.17*16^  b<.-).o+*  id3   (19) 

These  values  are  obviously  much  less  than  one,  but  the  factor  of  l/cos^5 
makes  the  coefficient  of  the  third  term  much  larger  near  the  pole.  At 
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80  N,  the  most  northerly  point  at  which  (15)  and  (16)  are  applied, 
cos   =  .174  so  b~/cos   =  0.047  still  well  within  the  region  of 
stability,  Unfortunately,  a  longer  time-step  drastically  reduces  the 
quality  of  the  results  and,  although  no  theoretical  investigation  is 
attempted,  this  may  be  caused  by  either  mass  imbalances  in  the  wind 
data  or  by  the  variability  of  the  space  increment  resulting  from  the 
cos^  effect. 

C.  WIND  DATA 

The  wind  data  used  by  Mahlman  (1967,  1969)  are  employed  in  the 
advective  terms  of  (11).  The  data,  which  cover  the  41  day  period  from 
15  November  until  25  December  1958,  consist  of  the  three  wind  components 
and  the  fields  of  temperature  and  geopotential  height  tabulated  at  50 
and  100  mb  for  the  region  from  40  to  80  N.  Originally  the  two  hori- 
zontal wind  components  as  well  as  the  temperature  and  height  fields 
were  extracted  directly  from  the  Weather  Bureau  (1963)  stratospheric 
daily  chart  series  for  the  IGY,  while  them's  were  computed  from  the 
other  fields  using  the  thermodynamic  equation  with  an  assumed  diabatic 
heating  rate  (Mahlman,  1967). 

During  the  early  portion  of  the  period  of  data  coverage,  the  polar 
night  vortex  is  intensifying  with  the  circulation  dominated  by  a  cold 
low  over  northern  Asia.  Then  on  21  November  the  flow  undergoes  a  pulsa- 
tion or  "minor  breakdown",  becoming  more  meridional  in  character  and 
remaining  disturbed  until  11  December  when  the  vortex  re-stabilizes  with 
an  important  wave  number  two  component.  During  this  last  period  one 
trough  lies  over  Asia  and  the  other  over  North  America  with  the  two 
intervening  ridges  over  the  oceans. 
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With  the  exception  of  the  ten-day  period  from  5  December  until  15 
December,  the  mean  cell  in  these  data  has  rising  motion  over  the  pole 
and  descent  in  middle  latitudes  (Mahlman,  1969).  Even  during  the  time 
in  December  when  the  vertical  motion  reverses  at  the  pole,  descent  takes 
place  only  in  the  immediate  region  of  the  pole  with  rising  motion 
continuing  in  the  60  to  70  N  band.  This  means  that,  except  for  that 
brief  period  in  December,  the  mean  cell  oeprates  in  the  exact  opposite 
sense  to  that  required  by  the  Brewer-Dobson  theory. 

Unfortunately,  the  data  are  tabulated  every   24  hours,  on  only  two 
pressure  surfaces,  and  for  a  somewhat  coarser  horizontal  grid.  This 
means  that  the  data  have  to  be  extended  by  extrapolation  and  interpolation 
in  both  space  and  time-  First,  using  the  thermal  wind  equation,  the  hor- 
izontal winds  are  extrapolated  upward  to  12.5  mb  and  downward  to  150  mb, 
and  the&Vs  at  these  two  levels  are  assumed  to  be  zero.  This  is  done 
at  every   point  in  the  horizontal  that  carries  the  initial  data  fields 
{every   40°  of  longitude  at  80  N,  every   20°  at  70  N,  and  every   10°  of 
longitude  at  50,  50  and  40  N),  Values  for  every   10°  longitude  at  70  N 
80  N  are  generated  by  linear  interpolation  between  the  tabulated  values 
and  the  same  process  is  applied  again  to  produce  winds  at  odd  multiples 
of  five  degrees  of  latitude.  Finally,  a  cubic  Lagrange  interpolating 
polynomial  is  fitted  to  the  data  at  each  horizontal  grid  intersection  for 
the  four  levels  and  evaluated  at  the  intermediate  pressure  levels  to 
complete  the  wind  field  for  each  day   The  wind  for  each  time  step  is 
obtained  from  the  daily  wind  fields  by  simple  linear  interpolation  in 
time. 

The  wind  data  require  still  more  treatment  before  they  are  suitable 
for  use  in  the  integration  scheme-  Due  to  inaccuracies  in  observation 
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and  processing,  the  probable  errors  in  the  v  wind  component  are  at  least 
as  large  as  the  real  mean  cell.  Furthermore,  the  omegas  generated  by 
the  thermodynamic  method  are  not  necessarily  consistent  with  the  diver- 
gences of  the  horizontal  wind.  Mass  imbalances  in  the  large-scale  eddies, 
while  undesirable,  tend  to  cancel  each  other  when  integrated  over  the 
entire  region,  but  a  mean  cell  with  spurious  mass  sources  may  seriously 
distort  the  results  of  the  computation.  It  was  decided  to  remedy  the 
situation  by  removing  the  natural  mean  cell  and  replacing  it  with  an 
artificial  time  independent  one  which  exactly  satisfied  the  zonally- 
averaged  continuity  equation.  The  following  transformation  accomplishes 

this  end: 

U"A      -    ut 


"A  \J<*  Y/4      \/* 


to'j,*  *    *>*&-*>$  +  f2-l/$     .  (20) 

Here  the  double  primed  quantities  are  the  new  values,  the  unprimed  the 
old,  the  tilded  the  zonal  mean  old  values,  and  the  capitals  the  contrived 
mean.  This  new  mean  is  usually  run  with  ascent  at  the  poles;  although 
some  experiments  were  run  with  the  cell  reversed  for  comparison.  Typical 
magnitudes  for  the  new  mean  were  on  the  order  of  a  knot  for  V  and  .5  mb 
day"  for -CI.  The  new  cell  was  constant  in  time  and  no  attempt  was  made 
to  introduce  time  variations  resembling  those  of  the  natural  mean  cell. 

While  use  of  an  artificial  mean  cell  seems  to  introduce  the  very   element 
of  arbitrariness  that  was  avoided  through  employing  observed  winds,  it 
should  be  pointed  out  that  this  approach  both  permits  experimentation  with 
mean  cells  of  differing  character  and  allows  use  of  mean  circulations  that 
are  more  or  less  consistent  with  heat  and  momentum  considerations. 
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III.   EXPERIMENTS 

A.  OZONE  SIMULATION 

1 .  Boundary  and  Initial  Conditions 

The  behavior  of  stratospheric  ozone  is  simulated  by  the  evolution 
of  an  inert  tracer  from  a  zonally  symmetric  initial  state  roughly  cor- 
responding to  the  observed  zonal  mean  ozone  distribution  (Hering,  1960). 
Throughout  the  integration  the  observed  initial  values  at  12.5  and  150 
mb  are  retained  at  the  top  and  bottom  boundary  points.  This  condition 
can  be  justified  since  it  is  assumed  that  the  value  at  12,5  mb  is  control- 
led only  by  photochemistry  and  is  time-independent  except  for  diurnal  and 
seasonal  variation  which  are  not  reproduced.  Similarly  the  mixing-ratio 
at  150  mb  is  assumed  to  represent  a  constant  equilibrium  between  down- 
ward transport  and  destruction  in  the  troposphere. 

Originally  it  was  the  intention  to  parameterize  advection  into 
the  region  from  the  tropics  by  retaining  the  observed  zonal  mean  values 
at  the  southern  boundary  points,  but  the  final  version  of  the  model  is 
permitted  to  generate  its  own  southern  boundary  condition  inasmuch  as 
the  points  at  35  N  are  all  set  to  the  generated  zonal  mean  of  those  at 
40  N. 

Near  the  pole,  the  cosine  of  the  latitude  becomes  small  causing 
an  apparent  singularity  in  the  prognostic  equation.  This  difficulty  can 
be  avoided  by  applying  a  cartesian  version  of  the  prognostic  equation  at 
the  pole,  but  a  simplier  solution  is  suggested  by  the  treatment  of  the 
southern  boundary,  so  the  values  at  85  N  are  set  to  the  zonal  mean  of 
those  at  80  N. 
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These  boundary  conditions  imply  that  the  system  is  not  closed 
since  boundary  fluxes  can  occur  through  all  surfaces  of  the  region  of 
integration.  A  horizontal  flux  through  the  southern  boundary  is  physically 
reasonable  since  advection  from  tropical  regions  is  an  important  source 
of  mid  and  high-latitude  ozone-  Similarly,  the  flux  over  the  pole  is 
to  be  expected;  although  the  boundary  scheme  used  does  not  properly 
simulate  it,  and  a  spurious  over-the-pole  transfer  arises  which  distorts 
the  final  results.  The  upper  boundary  condition  allows  ozone  to  be 
transported  downward  into  the  region  of  integration  but  does  not  permit 
it  to  leave  through  the  top.  The  values  at  12.5  mb  are  invariably  higher 
than  those  in  the  interior  of  the  region.  Since  the  vertical  motion  is 
small  near  the  upper  boundary,  diffusion  predominates  allowing  the  tracer 
to  flow  inward  but  never  outward.  A  similar  set  of  conditions  arises  at 
the  lower  boundary;  except  here  the  boundary  values  are  always  lower  than 
those  in  the  interior  so  only  outward  fluxes  can  occur.  Thus,  although 
no  source  of  sink  terms  are  explicitly  included  in  the  model,  the  boundary 
conditions  crudely  approximate  the  behavior  of  photochemical  sources  in 
the  upper  stratosphere  with  a  sink  at  the  tropopause. 
2.  Results 

Figures  1  through  5  are  a  series  of  charts  portraying  the  ozone 
mixing  ratio  and  geopotential  height  on  the  50  and  100  mb  levels  at 
intervals  of  5  days  during  the  period  of  integration.  Comparison  between 
the  50  mb  and  100  mb  pressure  levels  indicate  that  the  degree  of  continuity 
of  the  ozone  field  between  these  two  levels  is  comparable  to  that  observed 
in  the  geopotential  height  field  so  that  representation  of  more  than  one 
horizontal  section  of  the  ozone  field  is  superfluous,  and  only  the  100  mb 
surfaces  are  presented  after  25  November. 
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FIG.  7.  Simulated  zonal-mean  ozone  cross  section  in  jlgm-qm 
A:  Initial  cross  section;  B:  cross  section  for  20 
December  with  mean  ascent  at  the  pole;  C:  for  the 
same  date  with  a  zero  mean  cell;  and  D:  for  a  mean 
cell  with  descent  at  the  pole. 
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These  charts  show  a  number  of  qualitative  features  that  are 
reassuringly  similar  to  the  characteristics  of  the  real  atmospheric 
ozone  field.  Dobson's  (1924)  observations  of  ozone  maxima  in  the 
troughs  are  reproduced.  Similarly,  Bozkov's  (1968)  observations  of 
sinusoidal  isopleths  of  concentration  with  ozone  maxima  in  the  standing 
troughs  over  continents  and  minima  over  the  intervening  seas  appear  in 
these  fields. 

At  50  mb  and  60  N  the  normalized  correlation  coefficient  between 
cd  and  m,  computed  over  the  entire  period  of  integration,  is  0.22.  This 
means  that  the  high  mixing  ratios  in  the  lows  can  be  explained  in  terms 
of  advection  from  above  just  as  Reed  (1950)  and  Normand  (1953)  postulated. 

On  the  50  mb  surface  the  ozone  mixing  ratio  values  initially 
decrease  from  .053^gm/gm  at  35  N  to  0.028/^gm/gm  at  85  N  while  at  100 
mb  the  initial  values  increase  from  0.010  to  .031  in  the  same  interval. 
As  the  integration  proceeds  the  transport  processes  reverse  the  50  mb 
gradient  by  18  November  and  establish  a  mid-latitude  maximum  at  75  N 
by  25  November.  This  maximum  tends  to  drift  southward  appearing  most 
frequently  at  about  65  N  during  the  latter  portion  of  the  experiment.  It 
is  quite  persistent  in  this  latitude,  but  often  it  disappears  for  several 
days,  being  absent  on  a  total  of  19  of  the  41  simulated  days.  A  similar 
maximum  appears  in  the  100  mb  profile  on  2  December  but  it  is  less  pro- 
nounced and  not  nearly  so  persistent  in  time  (being  present  on  10  of  41 
days).  At  both  altitudes  the  maximum  zonal  mean  gradient  occurs  south 
of  the  latitude  of  maximum  mixing  ratio,  with  the  field  being  more  nearly 
flat  over  the  pole.  These  results  are  consistant  with  Bozkov's  work,  and 
since  they  arise  during  the  integration  and  do  not  represent  more  persist- 
ence of  the  initial  state,  they  enhance  confidence  in  the  model. 
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When  time  series  of  the  zonal  mean  mixing  ratio  and  of  the 
various  processes  modifying  it  are  considered  (see  FIG.  8),  it  is  found 
that  the  horizontal  eddies  and  horizontal  diffusion  are  very   highly  cor- 
related with  the  fluctuations  of  the  mixing  ratio.  The  other  processes, 
the  vertical  eddies,  the  vertical  and  horizontal  mean  cells  and  the 
vertical  diffusion  are  all  about  an  order  of  magnitude  smaller  than 
these  first  two  effects  and  did  not  produce  noticable,  short-term  changes 
in  the  mixing  ratio. 

The  direction  of  horizontal  eddy  transport  changes  during  the 
integration.  At  the  beginning  it  is  poleward,  but  during  the  pulsation 
of  the  polar  night  vortex  it  changes  sign  and  remains  negative  for  22 
days  until  the  vortex  begins  to  re-intensify.  When  the  normalized  mixing- 
ratio  and  meridional  velocity  correlation  coefficient  is  computed  for  the 
entire  period,  its  value  is  +  .023  at  50  mb  and  60  N.  But  when  the  data 
are  stratified  into  poleward  and  equatorward  transport  periods,  the  cor- 
relation during  the  former  is  +  .317  and  the  latter  -.237.  Since  the 
omega-mixing-ratio  correlation  is  +  0.22,  this  indicates  that  the  dif- 
ference in  the  effectiveness  between  the  vertical  and  horizontal  eddies 
is  primarily  caused  by  the  smaller  ratio  of  characteristic  velocity  to 
space  increment  for  the  former. 

Because  the  mean  cell  is  constant  in  time  and  the  zonal  mean 
gradient  is  nearly  so,  the  sign  and  magnitude  of  the  instantaneous  change 
resulting  from  the  mean  cell  is  quite  persistent.  This  means  that,  while 
the  short-term  effect  is  small,  the  mean  cell  produces  a  marked  change  in 
the  field  when  integrated  over  41  days. 

This  effect  is  well  illustrated  by  the  differences  among  the  20 
December  100  mb  fields  in  FIG.  5  and  FIG.  6.  Figure  6  shows  the  20 
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December  mixing-ratio  field  resulting  from  a  zero  mean  cell  as  well 
as  that  resulting  from  a  reversed  mean  cell.  Comparison  among  these 
three  charts  indicates  that  in  all  cases  the  features  are  very  similar 
in  shape  and  location  but  the  values  differ  as  though  the  fields  were 
translated  by  addition  of  a  constant.  This  is  \/ery   nearly  what  happened. 
Figure  7  shows  the  zonal  mean  cross  sections  for  each  of  these  three 
cases  on  25  December  as  well  as  the  initial  zonal  mean.  Each  different 
mean  cell  changes  the  inclination  of  the  mixing  ratio  isopleths,  tending 
to  increase  their  slopes  when  sinking  occurs  at  the  pole  and  decrease 
them  when  rising  takes  place  at  the  pole.  This  effect  is  most  important 
in  the  northern  portions  of  the  region  where  it  moves  the  lines  in  the 
vertical  while  toward  south  they  remain  at  about  the  same  altitude  in 
all  cases.  It  should  be  pointed  out  that,  while  this  effect  becomes 
important  only  after  long  integration  times,  it  will  eventually  attain 
an  equilibrium  with  the  other  processes  and  not  accumulate  indefinitely. 
Apparently  the  model  seeks  to  attain  just  such  an  equilibrium 

during  this  simulation.  At  50  mb  and  60  N  the  mixing  ratio  rises 

-1  i  \ 

asymptotically  to  a  value  of  about  60  juqm/gm       (See  FIG=  8).  The 

horizontal  eddies,  the  vertical  diffusion,  and  the  horizontal  mean  all 
appear  to  contain  decaying  transients  that  largely  damp  out  by  10  December. 
Thus,  apparently  about  25  days  are  required  to  recover  from  the  dis- 
equilibrium introduced  when  the  model  is  initiated  from  a  zonally  sym- 
metric initial  state, 

The  contribution  of  both  components  of  the  sub-grid-scale  diffusion 
is  disturbingly  large.  This  phenomenon  is  dependent  upon  a  completely 
arbitrary  sub-grid-scale  parameterization  scheme  and  probably  is  much 
stronger  than  any  real  atmospheric  diffusion.   It  is  unfortunate  that 
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the  diffusivities  can  not  be  reduced  without  impairing  the  model's 
computational  properties.  The  un-naturally  strong  diffusion  may  mask 
important  advective  processes  and  spuriously  shorten  the  equilibrium 
time  as  well  as  displace  the  equilibrium  from  its  true  value. 

B.  RADIOACTIVE  DEBRIS  SIMULATION 

1 .  Boundary  and  Initial  Conditions 

Clouds  of  radioactive  debris  are  simulated  by  permitting  the 
distribution  of  tracer  to  evolve  from  an  initial  three  dimensional 
Gaussian  distribution  centered  at  a  point  within  the  region  of  integra- 
tion. It  would  be  more  satisfying  to  use  an  initial  state  in  which  all 
the  tracer  is  concentrated  at  a  single  grid  point,  but  in  such  experi- 
ments an  undesirable  computational  mode  is  excited  by  the  steep  gradients 
so  that  the  results  are  meaningless. 

The  boundary  conditions  imposed  upon  the  evolving  cloud  of  debris 
represent  a  serious  problem.  For  the  previous  experiment  with  zonally 
symmetric  ozone  fields  as  an  initial  state,  the  generated  zonal  mean  of 
points  along  the  next  interior  latitude  circle  represents  a  reasonable 
boundary  value.  Because  the  point  source  problem  is  asymmetrical,  con- 
straining the  gradient  at  the  southern  boundary  to  zero  seems  to  constitute 
a  reasonable  boundary  condition,  but  computational  difficulties  arise  in 
this  case  when  the  wind  normal  to  the  boundary  is  strong.  The  final 
solution  to  this  problem  lies  in  a  compromise  between  the  two  approaches. 
Each  boundary  point  at  35  N  or  85  N  is  set  to  a  weighted  average  of  the 
zonal  mean  at  the  next  interior  parallel  and  the  value  at  the  adjacent 
point  along  that  parallel.  The  weighting  factors  are  0.7  times  the  mean 
and  0.3  times  the  adjacent  point  at  the  southern  boundary^ and  equal  weights 
are  used  at  the  northern  boundary. 
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FIG.  11.  Zonal-mean  simulated  radioactive  debris  in  arbitrary 
units  on  A:  18  November;  B:  20  November;  C:  22 
November;  D:  24  November;  and  D:  26  November. 
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At  the  top  and  bottom  of  the  region  there  is  little  inflow  or 
outflow  because  the  vertical  velocities  are  small.  Furthermore,  these 
regions  are  fairly  remote  from  the  centers  of  concentration.  This  means 
that  zero  gradient  at  12.5  and  at  150  mb  constitutes  a  physically  reason- 
able and  numerically  stable  boundary  condition  there, 
2.  Results 

The  evolution  of  a  cloud  injected  at  latitude  60  N  longitude 
140  W  with  initial  standard  deviation  10°  of  latitude,  20°  of  longitude 
and  25  mb  of  pressure  was  simulated  for  eleven  days  (See  FIGS.  9,  10). 
This  initial  position  lies  in  the  almost  purely  zonal  flow  in  the  ridge 
over  the  Northern  Pacific  so  that  the  main  mass  of  debris  is  carried 
eastward  along  the  height  contours  while  gradually  spreading  and  dis- 
persing. By  23  November  the  cloud  has  migrated  about  90°  of  longitude 
and  lies  in  the  eastern  portion  of  the  Atlantic  ridge.  During  this  period 
the  mixing  ratio  at  the  last  closed  contour  decreases  from  80,  its  initial 
value,  to  20  (arbitrary  units). 

On  17  November  part  of  the  cloud  is  carried  over  the  pole  and 
then  southeastward  around  the  Siberian  low.  Between  19  and  21  November 
the  mass  rounds  the  trough  and  begins  to  move  northward  so  that  by  23 
November  it  has  circumnavigated  the  pole  at  high  latitudes  and  re-merged 
with  the  main  cloud. 

In  the  zonal  mean,  the  center  of  the  initial  distribution  lay  at 
60  N  and  a  pressure  altitude  of  75  mb  (See  FIG.  11).  As  the  cloud  evolves 
the  center  of  maximum  mixing  ratio  gains  altitude  and  moves  toward  the 
south,  at  the  same  time  spreading  laterally  in  such  a  manner  that  its 
axis  slopes  downward  toward  the  pole.  By  23  November  the  zonal  mean 
bears  a  striking  resemblance  to  both  the  results  of  parameterized  studies 
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and  to  observations  (Guchksen  et  al . ,  1969),  showing  that  this  model  can, 
through  advection  and  diffusion  alone,  quickly  shape  any  reasonable  dis- 
tribution of  inert  tracer  into  a  fair  reproduction  of  observed  mixing 
ratio  fields.  In  fact,  this  process  took  place  much  more  quickly  than 
it  would  have  in  the  real  atmosphere  because  of  the  large  standard 
deviation  in  the  initial  distribution. 

No  detailed  analysis  is  made  of  the  principal  agencies  responsible 
for  the  tracer's  evolution,  but  the  horizontal  eddies  and  diffusion  com- 
bined are  about  an  order  of  magnitude  stronger  than  the  sum  of  the  other 
four  agencies.  However,  as  in  the  ozone  experiments,  the  sign  of  the 
tendencies  produced  by  the  mean  cell  terms  remain  consistent  in  time  so 
that  over  periods  longer  than  eleven  days  their  effect  becomes  significant 

IV,  CONCLUSIONS 

This  model  successfully  simulates  many  observed  qualitative  features 
of  stratospheric  transports  of  both  ozone  and  radioactive  debris.  Apart 
from  some  numerical  difficulties  which  result  from  mass  imbalance  in  the 
wind  field  and  from  inadequate  treatment  of  the  singularity  at  the  pole, 
no  obstacles  exist  to  prevent  the  use  of  this  model  to  predict  the  evolu- 
tion of  stratospheric  tracers. 

If  the  erroneous  exaggeration  of  the  parameterized  diffusion  may  be 
neglected,  the  instantaneous  derivative  is  primarily  determined  by  the 
horizontal  eddies,  Although  the  effect  of  the  vertical  eddies  is  about 
an  order  of  magnitude  weaker  than  that  of  the  horizontal  eddies,  it  can- 
not be  neglected  in  the  zonal  meanc  The  two  components  of  the  mean  cell 
are  still  weaker  than  the  vertical  eddies.  However,  over  long  periods  of 
time,  they  may  assume  greater  importance  because  they  do  not  change  sign 
or  fluctuate  greatly  in  magnitude. 
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